% main_multisector_planner
% Jones, Philippon, Venkateswaran (2021)
%

clear 
clc

%%

warning('Off','all') ;

path(pathdef) ;

%% Parameters

e0exog   = 0 ;
e0shift  = 0 ;
Imax     = 0.2;              %% targets for mortality
BaseMort = 0.0025;           
PeakMort = 3*BaseMort ;
deterministic_vac = 1 ;
stoch_prob = 155/156 ;

% Setup default sector-specific parameter values here
num_sectors = 4 ;
for j=1:num_sectors
	eval(['ec', int2str(j), ' = e0exog + (1-e0exog)*1/4 ;']) ;
	eval(['el', int2str(j), ' = e0exog + (1-e0exog)*1/4 ;']) ;
	eval(['chibar', int2str(j), ' = 0.5 ;']) ;
	eval(['Deltachi', int2str(j), ' = 0.34 ; ']) ;
	eval(['z', int2str(j), ' = 1 ; ']) ;
end

% Change sector-specific parameters here
chibar2 = .5*(.19/.45) ; % Easy to WFH
chibar3 = .5*(.76/.45) ; % High inability to WFH
chibar4 = .5*(.76/.45) ; % High inability to WFH

el2 = el1*(.19/.45) ;
el3 = el1*(.76/.45) ;
el4 = el1*(.76/.45) ;

ec2 = ec1*(.19/.45) ;
ec3 = ec1*(.76/.45) ;
ec4 = ec1*(.76/.45) ;

Deltachi2 = Deltachi1*(.19/.45) ; % Easy to WFH
Deltachi3 = Deltachi1*(.76/.45) ; % High inability to WFH
Deltachi4 = Deltachi1*(.76/.45) ; % High inability to WFH

beta       = (0.97)^(1/52)*stoch_prob ; 
ee_T_vac   = zeros(1,1000); 
if deterministic_vac; beta = (0.97)^(1/52); T_vac = 1.5*52 ; ee_T_vac(T_vac:end)=1;  end
sigma      = 4 ;
e0         = e0exog        + (1-e0exog)*1/2 ;
R0         = 2 ;
kappa      = 0.15;                                    % Fraction severely sick
rho        = 0.35;                                    % Recovery rate of infected 
deltabar   = BaseMort/(1-BaseMort)*rho/kappa;         % baseline death rate 
phi        = log(1-deltabar+PeakMort/(1-PeakMort)*rho/kappa)/(Imax*kappa);
N          = 1 ;
eta        = 2 ;
rhom       = 1  ;
us         = 0.5 ;
ud         = 2.5 ;
ee_i0_init = (1/1000)*N; % Initial shock of population
psim       = 0.99 ;
deltascale = 1 ;

model_name = 'infections_multisector_planner' ;
set_params_multisector_planner

run_dynare
 
%% For Figures

table = [ ...
    (1/1)*(L1(1:50))./(L4(1:50)+L1(1:50)+L2(1:50)+L3(1:50)), ...
    (1/1)*L2(1:50)./(L4(1:50)+L1(1:50)+L2(1:50)+L3(1:50)), ...
    (1/1)*(L3(1:50)+L4(1:50))./(L4(1:50)+L1(1:50)+L2(1:50)+L3(1:50)), ...
    (1/4)*(1-M1(1:50)).*(1-M1(1:50)).*el1.*L1(1:50).*L1(1:50)./e_by_l(1:50), ...    
    (1/4)*(1-M2(1:50)).*(1-M2(1:50)).*el2.*L2(1:50).*L2(1:50)./e_by_l(1:50), ...
    (2/4)*(1-M3(1:50)).*(1-M3(1:50)).*el3.*L3(1:50).*L3(1:50)./e_by_l(1:50)]' ;

table(:,13)

table2 = table(4:6,:)./table(1:3,:)-1 ;
table2(:,13) 

table3 = [M1' ; M2'; M3'] ;
4*table3(:,13)/(2*table3(1,13)+table3(2,13)+table3(3,13))-1

% Data. Avg Dingel-Neiman, Avg Employment, Avg Infections
% Columns are medium cost, low cost, high cost
data = [...
    0.40	0.25	0.22 
    0.76	0.30	0.13
    0.19	0.47	0.64];

tmp = cumsum(table([4,5,6],:)');

cumulative_infection = tmp./repmat(sum(tmp,2),[1, length(tmp(1,:))]) ; 
cumulative_infection = cumulative_infection(27,:) ; % After 26 weeks

%%

table = [ ...
    (1/1)*(C1(1:50)), ...
    (1/1)*C2(1:50), ...
    (1/2)*(C3(1:50)+C4(1:50))]' ;

cumulative_consumption = cumsum(table([1,2,3],:)');
cumulative_consumption = cumulative_consumption./...
    repmat((1:length(cumulative_consumption(:,1)))',[1, length(cumulative_consumption(1,:))]);
cumulative_consumption = cumulative_consumption(27,:) ;

table = [ ...
    (1/1)*(L1(1:50)), ...
    (1/1)*L2(1:50), ...
    (1/2)*(L3(1:50)+L4(1:50))]' ;

cumulative_labor = cumsum(table([1,2,3],:)');
cumulative_labor = cumulative_labor./...
    repmat((1:length(cumulative_labor(:,1)))',[1, length(cumulative_labor(1,:))]);
cumulative_labor = cumulative_labor(27,:) ;

[
    cumulative_infection
    cumulative_consumption
    cumulative_labor ]

[
    cumulative_infection - cumulative_infection(1)
    cumulative_consumption - cumulative_consumption(1)
    cumulative_labor - cumulative_labor(1)]

%%

load ./output/data_v_model.mat

data_v_model(:,4:6) = [
    cumulative_infection
    cumulative_consumption
    cumulative_consumption ];

data_v_model = data_v_model(:,[2 1 3 5 4 6]) ;

data_washington = [
0.944182150998049	1.01681588345397	0.967996593978898
0.883840115824324	0.930365061590146	0.930365061590146
0.888036697375743	0.919009671389634	0.919009671389634]';
data_v_model(2:3,1:3) = data_washington(2:3,:);

figure;

set(gcf, 'DefaultAxesLineWidth', 1);
set(gcf, 'DefaultLineLineWidth', 2);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',12);
set(gcf, 'Position', [1.6353e+03 712 560*(2) 420*(3/4)])

subplot(1,3,3) ;
toplot = reshape(data_v_model(3,:),[3 2]) ;
bar(toplot)
ylim([0.7 1.1]) ;
h=legend('Data','Model','location','northeast');
set(h,'Interpreter','latex')
xticklabels({'Low','Medium','High'}) ;
xlabel('Industry Group Infection Risk','Interpreter','latex')
title('GDP Ex-Tech, Avg at 26 Weeks','Interpreter','latex') ;

subplot(1,3,2) ;
toplot = reshape(data_v_model(2,:),[3 2]) ;
bar(toplot)
ylim([0.7 1.1]) ;
h=legend('Data','Model','location','northeast');
set(h,'Interpreter','latex')
xticklabels({'Low','Medium','High'}) ;
xlabel('Industry Group Infection Risk','Interpreter','latex')
title('GDP, Avg at 26 Weeks','Interpreter','latex') ;

subplot(1,3,1) ;
toplot = reshape(data_v_model(1,:),[3 2]) ; 
bar(toplot); hold on; 
plot([.75 1.25],[.25 .25], 'k','LineStyle',':')
plot([1.75 2.25],[.25 .25], 'k','LineStyle',':')
plot([2.75 3.25],[.5 .5], 'k','LineStyle',':')
xticklabels({'Low','Medium','High'}) ;
xlabel('Industry Group Infection Risk','Interpreter','latex')
title('Infections, \%','Interpreter','latex') ;
h=legend('Data','Model','SS Emp Share','location','northwest') ;
set(h,'Interpreter','latex') ;
ylim([0 0.7]) ;

print -depsc ./output/sector_bar.eps
print -dpng ./output/sector_bar.png
close
